Discrete Event Simulation in R

R packages for discrete event simulation modelling

Components of a simmer DES model

simmer models are built around the concept of a trajectory (i.e., an activity chain through which entities pass). In base simmer an action within a trajectory may be coded as below where a patient (the entity) is ‘seized’ by a nurse (the resource) to undertake an activity (e.g., a care process) for a given amount of time before releasing the entity, freeing up the resource.

trajectory("patients' path") %>%
  seize("nurse", 1) %>%
  timeout(function() rnorm(1, 15)) %>%
  release("nurse", 1)

We can however use the extension simmer.bricks to simplify this model which introduces the visit function which combines all of these steps into one function.

trajectory("patients' path") %>%
  visit("nurse", function() rnorm(1, 15))

Note that when part of the trajectory becomes more complex (for instance if there is a delayed release to allow for a resource to recover before it is used again) the base simmer solution may need to be reverted to, so it is important to be familiar with it.

A trajectory has a number of components that record activities that occur as a trajectory is followed.

Modelling a A&E department

To demonstrate a basic DES model, we can create a model of patients attending A&E. In the model a patient arrives and then the following trajectory is followed:

• The patient is booked in by a receptionist. The system has one receptionist and each booking in process takes around 5 minutes (modelled as a distribution rather than constant).

• The patient is then triaged by a nurse. The system has two nurses, and each triage process takes around 15 minutes (again modelled as a distribution). The nurse has to write up these findings and consequently is unable to see another patient for 5 minutes after the former patient has left.

• Finally, the patient is seen by a doctor to discuss ongoing care. The system has two doctors, and this process takes around 25 minutes.

We can then run this model by creating a generator which adds a new patient around every 5 minutes. We run the model for four hours.

set.seed(42)

# Create a simmer trajectory with the times taken for each action
# Note for nurse we cannot use visit as there is a delayed release of resource
patient <- trajectory("patients' path") %>%
  visit("receptionist", function() rnorm(1, 5)) %>%
  seize("nurse") %>%
  timeout(function() rnorm(1, 15)) %>%
  delayed_release("nurse", 5) %>% 
  visit("doctor", function() rnorm(1, 25))

# Create a simmer environment with the resources we have available
env <- simmer("A&E Department") %>%
  add_resource("receptionist", 1) %>%
  add_resource("nurse", 2) %>%
  add_resource("doctor", 2) %>%
  add_generator("patient", patient, function() rexp(1, 1/5))

# Run model for 4 hours (simmer does not have a time element. Our appointments are measured in minutes this is 4 * 60 minutes)
env %>% run(4 * 60)

## simmer environment: A&E Department | now: 240 | next: 242.578946591625
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 1(1) | queue status: 4(Inf) }
## { Resource: nurse | monitored: TRUE | server status: 1(1) | queue status: 17(Inf) }
## { Resource: doctor | monitored: TRUE | server status: 2(2) | queue status: 5(Inf) }
## { Source: patient | monitored: 1 | n_generated: 47 }

We can extract monitored data from the simulation using get_mon_resources. This allows us to see at different time points what each resource is doing.

resources <- get_mon_resources(env)
head(resources)

##       resource      time server queue capacity queue_size system limit
## 1 receptionist 0.9916841      1     0        1        Inf      1   Inf
## 2 receptionist 4.2961603      1     1        1        Inf      2   Inf
## 3 receptionist 6.3548125      1     0        1        Inf      1   Inf
## 4        nurse 6.3548125      1     0        2        Inf      1   Inf
## 5 receptionist 6.6620435      1     1        1        Inf      2   Inf
## 6 receptionist 8.8531660      1     2        1        Inf      3   Inf
##   replication
## 1           1
## 2           1
## 3           1
## 4           1
## 5           1
## 6           1

simmer.plot allows us to quickly visualise the results of the simulation. For instance, we can explore the utilisation of each resource.

plot(resources, metric = "utilization")

We can also see this over time. The blue line reports the number of staff busy, and the red line reports the number of patients queuing and waiting to be served. We could also show a line by adding “system” to the items argument to show the total number queuing and being served. This shows that for this simulation nurses are quickly overwhelmed (their capacity varies because of the admin work they have to do after every patient).

plot(resources, metric = "usage", items = c("queue", "server"))

It can sometimes to be easier to see this using the steps argument.

plot(resources, metric = "usage", items = c("queue", "server"), steps = T)

Finally, we can explore the waiting times (time in activity and total time in the system are also able to be plotted using “activity_time” and “flow_time” respectively.)

plot(get_mon_arrivals(env), metric = "waiting_time")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

We can also plot the trajectory itself which can be useful for communication.

plot(patient)

set.seed(42)

# Create a simmer environment with the resources we have available and run this for 4 hours. Replicate the model 100 times.
envs <- lapply(1:100, function(i) {
  simmer("A&E Department") %>%
    add_resource("receptionist", 1) %>%
    add_resource("nurse", 2) %>%
    add_resource("doctor", 2) %>%
    add_generator("patient", patient, function() rexp(1, 1/5)) %>%
    run(4 * 60)
})

plot(get_mon_resources(envs), metric = "usage", items = "queue")

plot(get_mon_arrivals(envs), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

set.seed(42)

# Create a simmer environment with the resources we have available and run this for 4 hours. Replicate the model 100 times.
envs <- lapply(1:100, function(i) {
  simmer("A&E Department") %>%
    add_resource("receptionist", 1) %>%
    add_resource("nurse", 4) %>%
    add_resource("doctor", 2) %>%
    add_generator("patient", patient, function() rexp(1, 1/5)) %>%
    run(4 * 60)
})

plot(get_mon_resources(envs), metric = "usage", items = c("server", "queue"))

plot(get_mon_arrivals(envs), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

set.seed(42)

# Create a simmer environment with the resources we have available and run this for 4 hours. Replicate the model 100 times.
envs <- lapply(1:100, function(i) {
  simmer("A&E Department") %>%
    add_resource("receptionist", 1) %>%
    add_resource("nurse", 3) %>%
    add_resource("doctor", 3) %>%
    add_generator("patient", patient, function() rexp(1, 1/5)) %>%
    run(4 * 60)
})

plot(get_mon_resources(envs), metric = "usage", items = c("server", "queue"))

plot(get_mon_arrivals(envs), metric = "waiting_time")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Optimisation

If we want to test how to maximise patient flow, we can use simmer.optim to iterate over a large number of models and assess what level of resource is optimal. Take the following:

• a maximization of the number of patients served before the 4-hour mark has been reached (**objective**)

• an employee budget of £2500 is available for the shift (**constraint**)

  • receptionist cost: £25 per hour

  • nurse cost: £50 per hour

  • doctor cost: £100 per hour

• a maximum average waiting time of 30 minutes is allowed (**constraint**)

We want to test having anywhere between 1 and 4 of each resource (receptionist, nurse, doctor). To do this manually would require 64 runs of the model but with simmer.optim we can do it in one run.

# remotes::install_github("r-simmer/simmer.optim")
library(simmer.optim)

# Create the model
model <- function(){
  patient <- trajectory("patients' path") %>%
    visit("receptionist", function() rnorm(1, 5)) %>%
    seize("nurse") %>%
    timeout(function() rnorm(1, 15)) %>%
    delayed_release("nurse", 5) %>% 
    visit("doctor", function() rnorm(1, 25))

  env <- simmer("A&E Department") %>%
    add_resource("receptionist", .opt("nr_receptionists")) %>%
    add_resource("nurse", .opt("nr_nurses")) %>%
    add_resource("doctor", .opt("nr_doctors")) %>%
    add_generator("patient", patient, function() rexp(1, 1/5))

  env
}

# Create the constraints
assert_within_budget <- function(envs, budget){
  runtime_hrs <- msr_runtime(envs) / 60
  number_receptionists <- msr_resource_capacity(envs, "receptionist")
  number_nurses <- msr_resource_capacity(envs, "nurse")
  number_doctors <- msr_resource_capacity(envs, "doctor")
  
  (number_receptionists * runtime_hrs * 25 +
    number_nurses * runtime_hrs * 50 +
    number_doctors * runtime_hrs * 100 ) <= budget
}

output <- simmer_optim(
  model = model,
  method = grid_optim,
  direction = "max",
  objective = msr_arrivals_finished,
  constraints = list(with_args(assert_within_budget, budget = 2500),
                     with_args(assert_avg_waiting_time_max, max_val = 30)),
  params = list(nr_receptionists = par_discrete(1:4),
                nr_nurses = par_discrete(1:4),
                nr_doctors = par_discrete(1:4)),
  control = optim_control(run_args = list(until = 4 * 60),
                          replications = 20))

output

## simmer.optim result 
## 
## method:                 grid_optim 
## objective value:        30.35 
## constraints satisfied:  TRUE 
## 
## params: 
## >  nr_receptionists :  1 
## >  nr_nurses :  4 
## >  nr_doctors :  4

We can plot the output of the simulation as before by referencing the envs object within the output list.

plot(get_mon_arrivals(output$envs), metric = "waiting_time")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot(get_mon_resources(output$envs), metric = "usage", items = c("server", "queue"))